Nonlinear photonic crystals near the supercollimation point 
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We uncover a strong coupling between nonlinearity and diffraction in a plio- 
tonic crystal at the supercoUimation point. We show this is modeled by a 
nonlinear diffraction term in a nonlinear schrodinger type equation, in which 
the properties of solitons are investigated. Linear stability analysis shows soli- 
tons are stable in an existence domain that obeys the Vakhitov-Kolokolov 
criterium. In addition, we investigate the influence of the nonlinear diffraction 
on soliton collision scenarios. 
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Photonic crystals 



PhC) are under active investigation due to the rich physics and all- 



optical signal control[l|, [2, y]. One o 
(SC) effect, that was experimentally 



' the striking features in PhCs is the supercollimation 



4, l5| and theoretically 



71, l8| examined. This effect 



originates from the possibility of obtaining fiat regions in the spatial dispersion relation or 
equifrequency contours of PhCs. At these particular points the phase propagation compo- 
nents along the direction of a beam are equal. In this way all components of the beam travel 



with the same phase velocity and it becomes nondiffr active. Recently, centimeter-scale SC 
was achieved in a large-area two-dimensional PhC 9|]. However, the study of nonlinear ef- 



fects around these SC points is largely unaddressed. Here we propose a more fundamental 



approach, as opposed to a more numerical one [10 1. 



In this letter we derive a phenomenological model that describes light beam propaga- 
tion in nonlinear PhCs around the SC point. In most optical systems, nonlinearity is a 
small perturbation, and hence added to the linear equation of motion. In contrast, in the 
current system, one uncovers a strong coupling between nonlinearity and diffraction. More 
specifically, a nonlinear diffraction term is introduced into the nonlinear Schrodinger (NLS) 
equation. This term can control the magnitude, and sometimes even the sign of the diffrac- 
tion. As a result, this is a unique system, and new physical phenomena emerge. For example, 
this modified equation gives rise to solitons up to an upper threshold for the propagation con- 
stant. Linear stability analysis demonstrates the solitons are stable in the existence domain, 
and obey the Vakhitov-Kolokolov criterium. In addition, we examine soliton interaction 
scenarios. 



In order to deduce the model equations, we consider a 3D PhC with the same periodicity 



in the transversal x and y directions (a scheme shown in Fig. [T]), and write the electric 
field of a beam in the PhC as E(x, y, z, t) = F{x, y, z)A{x, y, z) exp{ik'^z — iut), with k'^ the 
propagation constant of the central Bloch mode, uj the frequency of the beam, A{x, y, z) the 
slowly varying amplitude, and F(a;, y, z) the Bloch mode profile corresponding to uj and fc^. 
The propagation direction is along z. 

We study beams in the neighborhood of an SC point, so ~ ojsc with ojsc the SC 
frequency. The special dispersion relation around this point is shown in Fig. [H The sign 
and strength of diffraction depends strongly on the position of uj versus uJsc- Exactly at 
the frequency uJsci all components of the beam travel with the same propagation constant 
along z, noted by /cf", and there is no diffraction. For frequencies different than ujsc we 
approximate the equifrequency contours by a parabola. In this case diffraction can be 
positive or negative, depending on the position of u versus ujsc- Furthermore, the strength 
of diffraction increases as the difference between uj and uJsc increases. 



To include nonlinearity we apply first order perturbation theory Thus, the nonlinear 
interaction causes a small shift of the local dispersion relation [Fig. [1] , which is equivalent 
to shifting uJsc- To first order this shift is given by 

Aa; 1 /drn2(r)n(r)[(F ■ F)(F*-F*) + 2 |F|^] 



oJsc 4 /drn2(r) |F| 

^\A{x,y,z)\^ = K\A{x,y,z)\^ (1) 



with K a nonlinear coefficient calculated from the linear Bloch mode profile. 

The nonlinearity shifts usc to ujsc + so by Taylor expanding the dispersion relation 
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around the SC point, the modified dispersion relation becomes approximately 

/c^ - /cf^ + ai{uj - Use - 1^ \Afi^sc) 

= Piico -cosc-^i lAfcoscM + kl) (2) 

with kz the longitudinal, and k^., ky the transverse propagation vector components. The 
term with corresponds to the frequency dependence of the central /c^-component (thus at 
ky = 0). Similarly, the term with Pi describes the frequency change of the curvature (or 
the diffraction). Note that we neglect higher-order diffraction terms, such as fourth-order 
diffraction, as we are considering propagation of a broad beam with respect to the PhC 
period. We transform Eq. [2] into the space domain, and arrive at the following equation 
in dimensionless form: z|| + |aV^g — V^g + 7 {ql^ q = 0, where = ^ + 

with the transverse coordinates r] and ( scaled to the spatial characteristic width Wq and ^ 
is the longitudinal coordinate scaled to the free-space diffraction length = 2hWq/ A, for 



Gaussian-like beams 12|]. a = —2f3i{uj — uJsc)Ld/WQ indicates the linear diffraction strength, 
and its sign characterizes the type of linear diffraction. The novel nonlinear diffraction term 
is preceded by /? = 2j3iK,u!scc^Ld/WQ, with c the speed of light. The usual nonlinear term 
has 7 = —ainuscc^Ld. In the following, we employ Sn{I) > (thus n2 > 0) and Pi > 0. 

In this paper we will mainly show results concerning 2D PhCs, so the resulting equation 
takes the form 



.dq 1 d^q 1 ^ . .2 d^q . .2 , . 

Note that all the coefficients for this model can be deduced from rigorous numerical simu- 



lations 



91] . Recently, an equation similar to Eq. [3] was reported as the continuous approx- 



imation of the Salerno model[13|, [ij, llSj]. Here in contrast, we derive the model from a 
physical system, describing light beam propagation in nonlinear PhCs with SC. Eq. [3] con- 
serves the power U = / In a — (3 \qf dt]. The stationary solutions of Eq. Ohave the form 
liv^O = '?^('7) exp(i6(^), where a real function w{ri) and a real propagation constant b are 
found by iterative relaxation method. To analyze stability we examine perturbed solutions 
q = {w + u + iv)exp{ib^), where the real u{ri,^) and imaginary v{ri,^) perturbations can 
grow with complex rate 6 upon propagation. The eigenvalue problem linearized from Eq. [3] 
around w{rj) is solved numerically. 

One concludes that Eq. [3] allows soliton solutions [see e.g. Fig. [21(a)], but only in a fi- 
nite band of propagation constants. More specifically, there exists an upper threshold for 
the propagation constant (&co)) which depends on the strength of nonlinear diffraction (3 



[Fig. [2](c)]. Above this value peakon solutions 



13( 1 appear, which are unphysical given the as- 



sumptions of our physical model. To model what happens beyond that point in a real PhC, 
one would have to include more terms to Eq. [31 To emphasize, in our system it is the non- 
linear diffraction term that leads to the reduction of the semi-infinite band of propagation 
constants (as in a pure cubic NLS system) to a finite one, which is typical for soliton families 



in models with competing nonlinearities, such as the cubic-quintic NLS equationjlgj. It is 
noted that the power is a non-monotonic function of the propagation constant [FiglJI^b)]. 
Figlll^d) shows the dependence of the width of solitons on the amplitude (inset plot shows 
the dependence of maximum amplitude on the propagation constant). As one can see from 
this plot, the nonlinear diffraction has a significant effect at larger amplitudes. In addition, 
it should be pointed out that a lower power would be needed to generate a soliton in our 



system because the diffraction is so weak to start with. Thus our system may provide an 
experimentally favorable way to manipulate nonlinear waves. 

Linear stability analysis shows that solitons are stable in the whole domain of their ex- 
istence. An instability growth rate calculation is shown in the inset of Fig. M^c). It is 
noted that the Vakhitov-Kolokolov criterium applies to our system for fundamental soli- 
tons, namely solitons are stable when ^ > 0. To confirm the outcome of the linear stability 
analysis we perform direct numerical simulations of Eq. [31 We employ a split-step Fourier 
method, in combination with fourth-order Runge-Kutta to deal with the nonlinear diffrac- 
tion term. The input condition is q{i],^ = 0) = w{ri)[l + p{i])], with w{ri) the profile of 
the stationary soliton and p{ri) a random noise function with variance up to 10%. The 
simulations confirm the linear stability analysis. 

As another example where nonlinear diffraction affects fundamental phenomena, we re- 
port its influence on soliton collisions for different input conditions. We use as input a 
soliton of the pure NLS (/? = 0), and examine how changing f3 affects collision scenarios. 
In the case of two in-phase parallel solitons as input, the colliding solitons behave periodi- 
cally for a small nonlinear diffraction term [Fig. [3](a)], as in the case of the pure NLS. For 
larger /3 the colliding solitons merge into a single localized state, with breather-like features 
[Fig. [31(b)]. For solitons moving in opposite directions, the colliding solitons feature similar 
behavior. One example, presented in Fig. ^c), shows two solitons merging into a single lo- 
calized state. Such inelastic effects often appear when dealing with a perturbed NLS, again, 
here it is caused purely by nonlinear diffraction. For two out-of-phase solitons, the results 
show that the repulsive force between neighboring solitons is reduced with an increase of 
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nonlinear diffraction [Fig. [3](d)]. 

Finally, we investigated soliton properties in the (2+l)D model (two transversal dimen- 
sions and one propagation direction) but no stable solitons were found. The results show 
however that nonlinear diffraction modifies the critical power for collapse: positive (nega- 
tive) (3 reduces (increases) the critical power. This means the nonlinear diffraction can slow 
down the collapse. 

In most PhC structures the SC effect is relatively broadband. This means that the 
curvature changes slowly near the SC point. Therefore, to have a significant nonlinear 
diffraction effect one needs to design a PhC with a large curvature change, meaning a large 
(3i. At the SC frequency the nonlinear diffraction term would then be the main contribution 
that can interact with the normal nonlinear term. 
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Figure captions 

Figure 1. (a) A scheme of 3D PhC (b) Depiction of tlie linear dispersion relation in the 
proximity of a SC point, (c) Nonlinearity gives rise to an index change Sn, and is modeled 
by a shift of the dispersion relation, here shown for the particular input frequency cu — uJsc-- 

Figure 2. (a) Profile of a soliton for 6 = 1 with (3 — 0.3. (b) The dispersion diagram 
for 13 = 0.3. (c) Domain of existence for solitons in the (6, /3) plane, where the inset shows 
the real part of the perturbation growth rate versus propagation constant for /3 = 0.3. (d) 
The soliton width (FWHM) versus maximum amplitude for varying nonlinear diffraction, 
namely (3 = 0.1 (dotted line), 0.2 (dashed line) and 0.3 (solid line). The inset shows the 
dependence of maximum amplitude on the propagation constant, a — 1 and 7 = 1 for all 
cases. 

Figure 3. Collision scenarios between solitons with a — ^ — 1: Two in-phase parallel 
solitons for (a) P = 0.05 and (b) (3 = 0.15. (c) Two solitons moving in opposite directions 
with an angle {9 — 0.2), and P — 0.15. (d) Two out-of-phase parallel solitons for /3 = 0.15. 
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Fig. 1. (a) A scheme of 3D PhC. (b) Depiction of the hnear dispersion relation in the 
proximity of a SC point, (c) Nonhnearity gives rise to an index change and is modeled 
by a shift of the dispersion relation, here shown for the particular input frequency cu — uJsc- 
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Fig. 2. (a) Profile of a soliton for 6=1 with j3 = 0.3. (b) The dispersion diagram for 
P — 0.3. (c) Domain of existence for sohtons in the (6, /3) plane, where the inset shows 
the real part of the perturbation growth rate versus propagation constant for (3 — 0.3. (d) 
The soliton width (FWHM) versus maximum amplitude for varying nonlinear diffraction, 
namely P — 0.1 (dotted line), 0.2 (dashed line) and 0.3 (sohd hne). The inset shows the 
dependence of maximum amplitude on the propagation constant, a — 1 and 7 = 1 for all 
cases. 
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Fig. 3. Collision scenarios between solitons with a = 7 = 1: Two in-phase parallel solitons 
for (a) (3 — 0.05 and (b) (3 — 0.15. (c) Two solitons moving in opposite directions with an 
angle {9 — 0.2), and (3 — 0.15. (d) Two out-of-phase parallel sohtons for (3 — 0.15. 
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